#include"system_m.h"
#include"help.h"
#include<fstream>
#include<iostream>
#include<cmath>//abs,sqrt
#include<algorithm>//min,max

using std::ifstream;
using std::abs;
using std::min;
using std::max;


void System_M::setupOrbitals()
{
  vector<OrbitalType> OrbitalsTemp;
  vector<int> OrdersTemp;
  //split different m orbitals from System_Tz_SM
  for(int i=0;i<pSystem->Orbitals.size();i++)
    {
      int l=pSystem->Orbitals[i].l;
      int jj=pSystem->Orbitals[i].jj;
      int tz=pSystem->Orbitals[i].tz;
      DataType e=pSystem->Orbitals[i].e;
      for(int mm=-jj;mm<=jj;mm+=2)
	{
	  Orbitals.push_back(OrbitalType(l%2,mm,tz,e));
	  Orders.push_back(i);
	}
    }
  setupGroups();
}
  
void System_M::setupTwoBodyStates()
{
  TwoBodyStates.clear();
  for(int i=0;i<Orbitals.size();i++)
    for(int j=i+1;j<Orbitals.size();j++)
      {
	int Tz=(Orbitals[i].tz + Orbitals[j].tz)/2;
	int Par=(Orbitals[i].par + Orbitals[j].par)%2;
	int M=(Orbitals[i].mm + Orbitals[j].mm)/2;;
	TwoBodyStates.push_back(TwoBodyStateType(i,j,Tz,Par,M));
      }
  setupTwoBodyChannels();//noting that the twobody states are sorted after this calling.
}

const int  System_M::MatEledim=System_M::VSystem::MatEledim;//initialize of static variable
void System_M::setupTwoBodyMat()
{
  setupIndiceIn2BStates();
  int totalChannels=TwoBodyChannels.size();
  TwoBodyMat.resize(totalChannels);
  for(int channel=0;channel<totalChannels;channel++)
    {
      int num2BInChannel=TwoBodyChannels[channel].size();
      TwoBodyMat[channel].clear();
      TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
      for(int bra=0;bra<num2BInChannel;bra++)
	for(int ket=bra;ket<num2BInChannel;ket++)
	  {
	    int ab=TwoBodyChannels[channel][bra];
	    int cd=TwoBodyChannels[channel][ket];
	    int a=TwoBodyStates[ab].a;
	    int b=TwoBodyStates[ab].b;
	    int c=TwoBodyStates[cd].a;
	    int d=TwoBodyStates[cd].b;
	    TwoBodyMatEle MatEle(MatEledim);
	    get2BmatOnNeed(a,b,c,d,MatEle);
	    set2Bmat(channel,bra,ket,MatEle);
	  }
    }
}

void System_M::setupIndiceIn2BStates()
{
  int totalOrbitals=Orbitals.size();
  IndiceIn2BStates.resize((totalOrbitals*(totalOrbitals-1))/2);
  for(int i=0;i<TwoBodyStates.size();i++)
    {
      int a=TwoBodyStates[i].a;
      int b=TwoBodyStates[i].b;
      //a<b
      IndiceIn2BStates[mapab(a,b)]=i;
    }
}

void System_M::get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
{
  MatEle.setZero();
  if(a==b || c==d) return;//pauli principle
  int phase=1;
  //make sure a<b,c<d,(ab)<=(cd)
  if(a>b)
    {
      swap(a,b);
      phase*=-1;
    }
  if(c>d)
    {
      swap(c,d);
      phase*=-1;
    }
  int Mab=(Orbitals[a].mm + Orbitals[b].mm)/2;
  int Mcd=(Orbitals[c].mm + Orbitals[d].mm)/2;

  if( Mab!=Mcd ) return;
  int ab=IndiceIn2BStates[mapab(a,b)];
  int cd=IndiceIn2BStates[mapab(c,d)];
  if(ab==-1||cd==-1) return;
  int channelab=ChannelIndice[ab];
  int channelcd=ChannelIndice[cd];
  //channelab should = channelcd
  if(channelab!=channelcd) return;
  int bra=TwoBodyStateIndice[ab];
  int ket=TwoBodyStateIndice[cd];
  if(bra>ket)
    swap(bra,ket);
  get2Bmat(channelab,bra,ket,MatEle);
  MatEle*=phase;
}

void System_M::get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
{
  MatEle.resize(MatEledim);
  MatEle.setZero();

  if(a==b||c==d) return;//pauli
  int mma=Orbitals[a].mm;
  int mmb=Orbitals[b].mm;
  int Tzab=(Orbitals[a].tz+Orbitals[b].tz)/2;;
  int MMab=mma+mmb;
    
  int mmc=Orbitals[c].mm;
  int mmd=Orbitals[d].mm;    
  int Tzcd=(Orbitals[c].tz+Orbitals[d].tz)/2;
  int MMcd=mmc+mmd;

  if(Tzab!=Tzcd) return;
  if(MMab!=MMcd) return;
  if((Orbitals[a].par+Orbitals[b].par+Orbitals[c].par+Orbitals[d].par)%2) return;//parity check

  int MM=MMab;
  int ordera=Orders[a];
  int orderb=Orders[b];
  int orderc=Orders[c];
  int orderd=Orders[d];

  int jja=pSystem->Orbitals[ordera].jj;
  int jjb=pSystem->Orbitals[orderb].jj;
  int jjc=pSystem->Orbitals[orderc].jj;
  int jjd=pSystem->Orbitals[orderd].jj;

  int Jmin=max(abs(jja-jjb)/2,abs(jjc-jjd)/2);
  int Jmax=min((jja+jjb)/2,(jjc+jjd)/2);
  if(Jmin>Jmax) return;
  
  VSystem::TwoBodyMatEle Temp;
  for(int J=Jmin;J<=Jmax;J++)
    {
      pSystem->get2BmatAt(ordera,orderb,orderc,orderd,J,Temp);
      if(Temp.empty()) continue;
      if(ordera==orderb) Temp*=sqrt(2.);
      if(orderc==orderd) Temp*=sqrt(2.);
      MatEle+=Temp*cg(jja,jjb,J*2,mma,mmb,MM)*cg(jjc,jjd,J*2,mmc,mmd,MM);
    }
}


double System_M::jplus(int a,int b) const
{
  if(Orders[a]!=Orders[b])
    return 0;
  if( (Orbitals[a].mm-Orbitals[b].mm)/2 ==1 )
    {
      int jj=pSystem->Orbitals[Orders[b]].jj;
      int mm=Orbitals[b].mm;
      double j=0.5*jj;
      double m=0.5*mm;
      return sqrt( j*(j+1)-m*(m+1) );
    }
  else
    return 0;
}

double System_M::jminus(int a,int b) const
{
  if(Orders[a]!=Orders[b])
    return 0;
  if( (Orbitals[b].mm-Orbitals[a].mm)/2 == 1 )
    {
      int jj=pSystem->Orbitals[Orders[b]].jj;      
      int mm=Orbitals[b].mm;
      double j=0.5*jj;
      double m=0.5*mm;
      return sqrt( j*(j+1)-m*(m-1) );
    }
  else
    return 0;
}
  
void System_M::printOrbitals(int num) const
{
  if(num==-1) num=Orbitals.size();
  std::cout<<"num"<<"\t"<<"par"<<"\t"<<"mm"<<"\t"<<"tz"<<"\t"<<"e"<<"\t"<<endl;        
  for(int i=0;i<num;i++)
    {
      std::cout<<i<<"\t"<<Orbitals[i].par<<"\t"<<Orbitals[i].mm<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<std::endl;
    }
}


void System_M::printTwoBodyStates() const
{
  for(int i=0;i<TwoBodyStates.size();i++)
    {
      std::cout<<TwoBodyStates[i].a<<"\t"<<TwoBodyStates[i].b<<"\t"<<TwoBodyStates[i].Tz<<"\t"<<TwoBodyStates[i].Par<<"\t"<<TwoBodyStates[i].M<<std::endl;
    }
}
